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A  FurtpSn-TV' BUbroutiae  is  grVen  which 


A  FurtPan-TV' BUbroutiae  is  grVen  which  may  be  used  to  generate 
random  observations  of  a  continuous  random  variable  from  a  table  of  dis¬ 
crete  observations.  This  subroutine  employs  Akiina's  algorithm  for  cubic 
spline  interpolation. 

The  compatibility  of  problem  and  algorithm  is  such  that  this  simple 
routine  gives  very  good  results. 


1.  The  Algorithm.  An  Important  part  of  any 
simulation  study  is  the  sampling  of  a  ran¬ 
dom  variable  x  with  cummulative  distribu¬ 
tion  function  F (x) .  Often  F  is  specified 
only  by  a  table  of  (x,F(x)),  the  probability 
distribution  of  the  random  variable  y“F(x) 
is  the  uniform  distribution  U(0,1).  There 
are  satisfactory  pseudorandom  number  gen¬ 
erators  to  sample  from  U(0,1).  For  ex¬ 
ample  from  Knuth  [2] 


(1.1)  y„ 


^  31Al5926y^+  27182818(mod  2’' 


The  problem  of  samoling  from  F(x)  may  be 
looked  at,  then,  as  equivalent  to  finding 
an  inversion  formula 

(1.2)  x.F-l(y)  . 

To  approximate  F~^(y)  we  use  an  empiri¬ 
cally  obtained  table 

(1.3)  ((x^,F(x^)):x^  <X2  <  ...  <Xj^) 

plus  a  good  interpolation  routine.  He 
should  have  in  (1.3)  that  F(Xj^)'--0  and 
F(Xj^)~l  .  In  this  way  it  la  possible  to 
Lake  random  values  from  U(0,1)  and  obtain 
their  approximate  Inverses  as  random  samples 
from  "nearly  F(x)." 
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The  sampling  problem  now  reduces  to  one 
of  interpolation  subject  to  a  montonicity 
constraint  or  "isotone  interpolation." 

Neither  the  standard  polynomial  interpolation 
nor  the  cubic  spline  Interpolation  have  a 
chance  of  accomplishing  this  objective.  For 
this  reason  (and  others)  the  standard  approach 
in  the  literature  is  piecewise  linear  inter¬ 
polation.  Observe  that  piecewise  liner  inter¬ 
polation  maintains  monotonicity  but,  in 
general,  at  the  expense  of  accuracy. 

It  is  the  purpose  of  this  note  merely  to 
point  out  that  Akima's  (1}  quasl-Hermite 
piecewise  cubic  interpolation  should  be  very 
effective  on  this  problem;  since  it  was  de¬ 
signed  to  handle  similar  problems.  A  random 
number  generator  (RNG)  is  advocated  which 
uses  Aklma  interpolation.  The  user  reads  in 
a  tabic 

((*(i),F(X(i))):X(o)  <x^^^  <  ...  <*(fcfl)^' 

The  routine  then  computes  and  stores  the  co¬ 
efficients  of  x“F"^(y).  Calls  on  the  ran¬ 
dom  number  generator  will  first  obtain  a 
pseudorandom  number  v  from  U(0,1).  If 
y<F(*(i)).  the  random  number  given  is 

if  y>F(x^j^j),  the  random  number  given  is 
If  F(x^j  <y  <F(x.j^.  )  ,  the  random  number  given 


0 


CAUCHY  DISTRIBUTION 


IsF  firl .  A  listing  of  the  random  number  gen¬ 
erator  is  given  in  Appendix  11. 

Naturally,  the  quality  of  the  random 
number  generator  proposed  is  determined  by 
the  ability  of  F'l  to  approximate  a  var¬ 
iety  of  cumulative  distribution  functions. 

In  Appendix  1.  we  show  how  the  interpola¬ 
tion  errors  of  F'^(y)  compare  with  those 
of  linear  interoolation  in  the  case  of  the 
standard  normal,  Cauchy  and  double  exponen¬ 
tial  distributions  with  location  0  and 
scale  1.  Exact  values  of  F~^(y)  are 
assianed  given  for  y  between  .50  and  .98 
in  increments  of  .01  and  for  y  between 
.98  and  .995  in  increments  of  .005.  In¬ 
terpolated  values  of  and  linear  inter¬ 

polation  are  obtained  at  points  between  the 
mesh  values  and  the  relative  errors  com¬ 
puted  and  compared. 
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APPENDIX  I 
NORMAL  DISTRIBUTION 

F(x)-=  X-  E(AKIMAS)=  E  (LINEAR) 


0.50333 

0.00833 

0.52333 

0.05837 

0.5A333 

0.10858 

0.56333 

0.15910 

0.58333 

0.21005 

0.60333 

0.26157 

0.62333 

0.31381 

C. 64333 

0.36694 

0.66333 

0.42114 

0.68333 

0.47662 

0.70333 

0.53362 

0.72333 

0.59241 

0.74333 

0.65334 

0.76333 

0.71680 

0.78333 

0.78329 

0.80333 

0.85343 

0.82333 

0.92805 

0.84333 

1.00823 

0.86333 

1.09546 

0.88333 

1.19193 

0.90333 

1.30097 

’.•478I2 

0.94333 

1.58372 

0.96333 

1.79115 

0.98333 

2.12849 

0.99333 

2.47156 

0.00006 

0.00011 

0.00003 

0.00020 

0.00002 

0.00033 

0.00006 

0.00044 

0.00004 

0.00061 

0.00001 

0.00080 

0.00001 

0.00096 

0.00000 

0.00114 

0.00004 

0.00129 

0.00005 

0.00149 

0.00004 

0.00172 

0.00002 

0.00200 

0.00000 

0.00231 

0.00007 

0.00258 

0.00004 

0.00300 

0.00005 

0.00346 

0.00007 

0.00401 

0.00012 

0.00469 

0.00011 

0.00561 

0.00018 

0.00679 

0.00021 

0.00853 

0.00030 

0.01115 

0.00034 

0.01569 

0.00185 

0,02539 

0.00336 

0.02805 

0.02432 

0.07186 

F(x)  = 

x= 

E(AKIMAS)  = 

E(  LINEAR): 

0.50333 

0.01047 

0.00002 

0.00010 

0.52333 

0.07344 

0.00005 

0.00054 

0.54333 

0.13698 

0.00005 

0.00098 

0.56333 

0.20164 

0.00005 

0.00143 

0.58333 

0.26795 

0.00005 

0.00190 

0.60333 

0.33654 

0.00006 

0.00237 

0.62333 

0.40809 

0.00006 

0.00287 

0.64333 

0.48342 

0.00007 

0.00340 

0.66333 

0.56347 

0.00008 

0.00396 

0.68333 

0.64941 

0.00009 

0.00456 

0.70333 

0.74267 

0.00011 

0.00521 

0.72333 

0.84507 

0.00013 

0.00592 

0.74333 

0.95897 

0.00015 

0.00672 

0.76333 

1.08749 

0.00017 

0.00761 

0,78333 

1.23490 

0.00021 

0.00865 

0.80333 

1.40714 

0.00025 

0.00985 

0.82333 

1.61283 

0.00030 

0.01128 

0.84333 

1.86499 

0.00037 

0.01304 

0.86333 

2.18419 

0.00046 

0,01526 

0.88333 

2,60509 

0.00057 

0.01820 

0.90333 

3.19100 

0.00068 

0.02229 

0.92333 

4.07127 

0.00072 

0.02844 

0.94333 

5.55776 

0.00021 

0.03882 

0.96333 

8.64274 

0,00113 

0.06036 

0,98333 

19.08113 

0.02053 

0.06658 

0,99333 

47,73947 

0.08864 

0.16661 

DOUBLE  EXP 

DISTRIBUTION 

F(x)- 

X“ 

E(AKIMAS)- 

B(  LINEAR) 

0.50333 

2.09951 

0.00002 

0.00223 

0.52333 

2.22281 

0.00003 

0.00233 

0.54333 

2,35140 

0.00002 

0.00244 

0.56333 

2.48575 

0.00002 

0,00255 

0.58333 

2.62641 

0.00002 

0.00268 

0.60333 

2.77398 

0.00003 

0.00281 

0.62333 

2.92918 

0.00003 

0.00295 

0,64333 

3.09286 

0.00004 

0,00312 

0.66333 

3.26599 

0.00003 

0.00331 

0.68333 

3.44972 

0.00004 

0.00351 

0.70333 

3. S' 544 

0.00005 

0.00375 

0.72333 

3.85483 

0.00006 

0.00402 

0.74333 

4.07993 

0.00006 

0.00434 

0.76333 

4.32331 

0.00007 

0.00471 

0.78333 

4.58819 

0.00009 

0.00514 

0.80333 

4.87874 

0.00010 

0.00567 

0.82333 

5,20047 

0.00013 

0.00631 

0.84333 

5.56091 

0.00015 

0.00712 

0.86333 

5.97063 

0.00019 

0.00817 

0.88333 

6.44530 

0.00024 

0.00957 

0.90333 

7.00946 

0.00031 

0.011 

57 

0.92333 

7.70487 

0.00040 

0.01461 

0.94333 

8.61171 

0.00042 

0.01983 

0.96333 

9.91767 

0.00215 

0.03088 

0.98333 

12.28305 

0.00450 

0.03291 

0,99333 
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^^interpolation  pF_T_HI-  INVERSE  .NORML  OIST«IOUTfON 
COilMON  N  «P~t  f  t»0)  .P’2(  SO>«P3<  50  )  •PA  1 50  >  •  X I  5C  > 

<  N  IS  The  NUMBER  OK  INTERPOLATING  KNOTS*  N  SHOULD  OC  LESS  OR  EO  TO 
N  «  SO 
CALL  INIT 

C  M  IS  THE  NUHOER  OK  RANDOM  NUMSERS  TO  BE  CBNFPATEO  M  SMALLER  THAN 
H  nJiC 
'call  RNGIMI 

CALL  EXIT 
ENO 

susi^utine  iNijr 

c  ••••••••••••*•■•••• 

C  THIS  SUBROoriNE  CEnCMATES  TnE  PARAMETERS  OK  THE  CUSlCS  THAT 
C  INTERPOLATE  THE  DATA  KNOTS 

C  ARRAYS  Pt«P?*P3.P«  CONTAIN  THF  PARAMETERS 

COMMON  N.Pl  <50)  .l>?(50}  .PilSO)  .PA  1501  •AI50I 
100  KQHMAT  t2K12«6> 

200  KOUMAT  (2K«tOK12«5) 

c  #•••«««  REAO  DATA  •*«••••••••«•••••••••*••• 

READ  100*( (X< I) .PAI t) ) •! ») *N) 

c  «••«•••««•  READ  DATA  FINISHES  *••••••*••• 

N2  *  N-? 

)  )/<X<2l^X<l  )  ) 

S2*<  PA<  3  l  -PA  (  2)  I  /iVc  3)-X<Vr) 

S3a(P4( 3).P4(4» )/<X(3)-X<A)) 

OO  10  |s3.N2 

S4  «  <PA<I«2»  >PA(I«ll)/(X(tf2}  -  XClAin 
»1  «  A0S($2*S1> 

_WI_  A6S<  S4-S3) 
tF~'(R|Aw3>  20*30.20 
20  R3<tl  s  <«3*S2  ♦  Wl«S3l/tW)4W3) 

GO  TO  40 

SO  P3(l)  «  .SAfSEASSI 

40  St  «  S2 

S2  m  S3 
10  si  *«  S4 

00  00  t«3«N2 
Aa- X< (♦! >  -  X(  I  ) 

R|<ll  *  CPJTf#  ♦  P3<f*tf  -  2.04fP4fr*l/  -  P4I r Jl/Al/<A«A) 
R2<lt«(3«04(P4<t>|).t.4(t))  /A  -  2*04P3<I)  -  P3n4||)/A 
-SA  .SQNTlNyB 
RETURN 

End 

SUBROUTINE  RNG(M). 

COMMON  N«P1 ISO) ,021 90) .PSISO) vPAlSO) •XlSO) 

REALTS  RN<?0)»I34«I35 
too  format  (l6*9Ft2,6) 

c  444444  RANDOM  NUMBER  GENERATOR  «444aa44 
134  «  5772 I S66 
,135  Jf. 2. 04435 

b6”lo'|a| .M 

134  a  I  314IS926  4134  ♦  27I82B2B  > 

134  «  (DMaO<l34«l3S)) 

RNIII  a  134/135 
PRINT  t00«l«RNf|l 
IP  ,<RNf  fj  •LE«X(3))  RNfl)  a  XCSI 
10  IF  (RNrri*G€«X<N«2|  )  RNI  Ua  XfN-2) 

•  IC  ■  I 
NS  a  N-3 
00  20  X  a3«N3 
00  30  Jat*M 
V  a  RNIJ)  •  Xlll 
'if  iRNIJt-XCIll  40,S0*50 
OB  IF  CRNf^l  •  X(|4|||  «0»60*40 

00  VRN  a  Pt(l»4V4V4V  ♦  P2f||4V4V  ♦  P3in4V  ♦  P4|ll 

c  TRN  IS  THE  GENEOATEO  RANDOM  OBSERVATION 
PRINT  IOO*IC«aN< Jl« VRN 
IC  •  IC  ♦  t 


40 

CONTINUE 

IF  lIC.CT 

20 

COhIInuE 

«• 

continue 

• 

continue 

f^TURN 

ENO 


/ 

1,  This  ifork  wos  oupportod  In  port  by  NASA  contrsct  NAS“9**1277^, 
OKU  Crxxt  Ot’7.-042-283  and  ERPA  cor.trcct  E]^40-1)-5046,V^ 
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